!----------------------------------------------------------------------! 
!       read data for FVM                                              !
!----------------------------------------------------------------------!
        SUBROUTINE READER()
        USE BASIC_MOD
        IMPLICIT NONE
        INTEGER:: I,J
        REAL*8:: S1,S2

!-----read mesh file
        OPEN(3,FILE='210-80.x',FORM='UNFORMATTED')
        READ(3)
        READ(3) NI,NJ      
        ALLOCATE (X(-1:NI+2,-1:NJ+2),Y(-1:NI+2,-1:NJ+2)) 
        READ(3) ((X(I,J),I=1,NI),J=1,NJ),((Y(I,J),I=1,NI),J=1,NJ)  
        CLOSE(3) 
        XCELLS=NI-1
        YCELLS=NJ-1
	    DO J=1,NJ
	        DO I=1,NI
	            X(I,J)=X(I,J)*10.d0
	            Y(I,J)=Y(I,J)*10.d0
	        END DO 
        END DO 
        DX=X(2,1)-X(1,1)
        DY=Y(1,2)-Y(1,1)
        DO J=1,NJ
            X(0,J)=2.D0*X(1,J)-X(2,J)
            Y(0,J)=2.D0*Y(1,J)-Y(2,J)
            X(-1,J)=2.D0*X(1,J)-X(3,J)
            Y(-1,J)=2.D0*Y(1,J)-Y(3,J)
            X(NI+1,J)=2.D0*X(NI,J)-X(NI-1,J)
            Y(NI+1,J)=2.D0*Y(NI,J)-Y(NI-1,J)
            X(NI+2,J)=2.D0*X(NI,J)-X(NI-2,J)
            Y(NI+2,J)=2.D0*Y(NI,J)-Y(NI-2,J)
        END DO 
        DO I=1,NI
            X(I,0)=2.D0*X(I,1)-X(I,2)
            Y(I,0)=2.D0*Y(I,1)-Y(I,2)
            X(I,-1)=2.D0*X(I,1)-X(I,3)
            Y(I,-1)=2.D0*Y(I,1)-Y(I,3)
            X(I,NJ+1)=2.D0*X(I,NJ)-X(I,NJ-1)
            Y(I,NJ+1)=2.D0*Y(I,NJ)-Y(I,NJ-1)
            X(I,NJ+2)=2.D0*X(I,NJ)-X(I,NJ-2)
            Y(I,NJ+2)=2.D0*Y(I,NJ)-Y(I,NJ-2)
        END DO 
        X(0,0)=2.D0*X(1,1)-X(2,2)
        Y(0,0)=2.D0*Y(1,1)-Y(2,2)
        X(0,-1)=2.D0*X(1,1)-X(2,3)
        Y(0,-1)=2.D0*Y(1,1)-Y(2,3)
        X(-1,0)=2.D0*X(1,1)-X(3,2)
        Y(-1,0)=2.D0*Y(1,1)-Y(3,2)
        X(-1,-1)=2.D0*X(1,1)-X(3,3)
        Y(-1,-1)=2.D0*Y(1,1)-Y(3,3)
        X(NI+1,NJ+1)=2.D0*X(NI,NJ)-X(NI-1,NJ-1)
        Y(NI+1,NJ+1)=2.D0*Y(NI,NJ)-Y(NI-1,NJ-1)
        X(NI+1,NJ+2)=2.D0*X(NI,NJ)-X(NI-1,NJ-2)
        Y(NI+1,NJ+2)=2.D0*Y(NI,NJ)-Y(NI-1,NJ-2)
        X(NI+2,NJ+1)=2.D0*X(NI,NJ)-X(NI-2,NJ-1)
        Y(NI+2,NJ+1)=2.D0*Y(NI,NJ)-Y(NI-2,NJ-1)
        X(NI+2,NJ+2)=2.D0*X(NI,NJ)-X(NI-2,NJ-2)
        Y(NI+2,NJ+2)=2.D0*Y(NI,NJ)-Y(NI-2,NJ-2)
        X(NI+1,0)=2.D0*X(NI,1)-X(NI-1,2)
        Y(NI+1,0)=2.D0*Y(NI,1)-Y(NI-1,2)
        X(NI+2,0)=2.D0*X(NI,1)-X(NI-2,2)
        Y(NI+2,0)=2.D0*Y(NI,1)-Y(NI-2,2)
        X(NI+1,-1)=2.D0*X(NI,1)-X(NI-1,3)
        Y(NI+1,-1)=2.D0*Y(NI,1)-Y(NI-1,3)
        X(NI+2,-1)=2.D0*X(NI,1)-X(NI-2,3)
        Y(NI+2,-1)=2.D0*Y(NI,1)-Y(NI-2,3)
        X(0,NJ+1)=2.D0*X(1,NJ)-X(2,NJ-1)
        Y(0,NJ+1)=2.D0*Y(1,NJ)-Y(2,NJ-1)
        X(-1,NJ+1)=2.D0*X(1,NJ)-X(3,NJ-1)
        Y(-1,NJ+1)=2.D0*Y(1,NJ)-Y(3,NJ-1)
        X(0,NJ+2)=2.D0*X(1,NJ)-X(2,NJ-2)
        Y(0,NJ+2)=2.D0*Y(1,NJ)-Y(2,NJ-2)
        X(-1,NJ+2)=2.D0*X(1,NJ)-X(3,NJ-2)
        Y(-1,NJ+2)=2.D0*Y(1,NJ)-Y(3,NJ-2)
!-----confirm the size of dynamic array
        ALLOCATE (XC(0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (YC(0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (VOL(0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (AREA_S(3,0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (AREA_E(3,0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (AREA_N(3,0:XCELLS+1,0:YCELLS+1))
        ALLOCATE (AREA_W(3,0:XCELLS+1,0:YCELLS+1))
!-----calculate the volume and the direction of interface of each cell
        DO I=0,XCELLS+1
            DO J=0,YCELLS+1
                S1=0.5D0*DABS((X(I+1,J)-X(I,J))*(Y(I,J+1)-Y(I+1,J))-(Y(I+1,J)-Y(I,J))*(X(I,J+1)-X(I+1,J)))
                S2=0.5D0*DABS((X(I+1,J+1)-X(I,J+1))*(Y(I,J+1)-Y(I+1,J)) -(Y(I+1,J+1)-Y(I,J+1))*(X(I,J+1)-X(I+1,J)))
                XC(I,J)=(S1*(X(I,J)+X(I+1,J)+X(I,J+1))/3+S2*(X(I+1,J+1)+X(I+1,J)+X(I,J+1))/3)/(S1+S2)
                YC(I,J)=(S1*(Y(I,J)+Y(I+1,J)+Y(I,J+1))/3+S2*(Y(I+1,J+1)+Y(I+1,J)+Y(I,J+1))/3)/(S1+S2)
                VOL(I,J)=S1+S2 
                AREA_S(1,I,J)=DSQRT((X(I+1,J)-X(I,J))**2+(Y(I+1,J)-Y(I,J))**2)
                AREA_S(2,I,J)=-(Y(I,J)-Y(I+1,J))/AREA_S(1,I,J)
                AREA_S(3,I,J)= (X(I,J)-X(I+1,J))/AREA_S(1,I,J)
                AREA_E(1,I,J)=DSQRT((X(I+1,J+1)-X(I+1,J))**2+(Y(I+1,J+1)-Y(I+1,J))**2)
                AREA_E(2,I,J)=-(Y(I+1,J)-Y(I+1,J+1))/AREA_E(1,I,J)
                AREA_E(3,I,J)= (X(I+1,J)-X(I+1,J+1))/AREA_E(1,I,J)
                AREA_N(1,I,J)=DSQRT((X(I+1,J+1)-X(I,J+1))**2+(Y(I+1,J+1)-Y(I,J+1))**2)
                AREA_N(2,I,J)=-(Y(I+1,J+1)-Y(I,J+1))/AREA_N(1,I,J)
                AREA_N(3,I,J)= (X(I+1,J+1)-X(I,J+1))/AREA_N(1,I,J)
                AREA_W(1,I,J)=DSQRT((X(I,J+1)-X(I,J))**2+(Y(I,J+1)-Y(I,J))**2)
                AREA_W(2,I,J)=-(Y(I,J+1)-Y(I,J))/AREA_W(1,I,J)
                AREA_W(3,I,J)= (X(I,J+1)-X(I,J))/AREA_W(1,I,J)
            END DO
        END DO
        DO J=1,YCELLS
            XC(0,J)=X(1,J+1)+X(1,J)-XC(1,J)
            YC(0,J)=Y(1,J+1)+Y(1,J)-YC(1,J)
            XC(XCELLS+1,J)=X(XCELLS+1,J+1)+X(XCELLS+1,J)-XC(XCELLS,J)
            YC(XCELLS+1,J)=Y(XCELLS+1,J+1)+Y(XCELLS+1,J)-YC(XCELLS,J)
        END DO
        DO I=1,XCELLS
            XC(I,0)=X(I+1,1)+X(I,1)-XC(I,1)
            YC(I,0)=Y(I+1,1)+Y(I,1)-YC(I,1)
            XC(I,YCELLS+1)=X(I+1,YCELLS+1)+X(I,YCELLS+1)-XC(I,YCELLS)
            YC(I,YCELLS+1)=Y(I+1,YCELLS+1)+Y(I,YCELLS+1)-YC(I,YCELLS)
        END DO
        XC(0,0)=2*X(1,1)-XC(1,1)
        YC(0,0)=2*Y(1,1)-YC(1,1)
        XC(XCELLS+1,0)=2*X(XCELLS+1,1)-XC(XCELLS,1)
        YC(XCELLS+1,0)=2*Y(XCELLS+1,1)-YC(XCELLS,1)
        XC(XCELLS+1,YCELLS+1)=2*X(XCELLS+1,YCELLS+1)-XC(XCELLS,YCELLS)
        YC(XCELLS+1,YCELLS+1)=2*Y(XCELLS+1,YCELLS+1)-YC(XCELLS,YCELLS)
        XC(0,YCELLS+1)=2*X(1,YCELLS+1)-XC(1,YCELLS)
        YC(0,YCELLS+1)=2*Y(1,YCELLS+1)-YC(1,YCELLS) 
!-----check mech 

        OPEN(2,FILE='check_grid.dat')
          WRITE(2,*) 'VARIABLES="X","Y"'
          WRITE(2,*) 'ZONE I=',NI+4,',J=',NJ+4
          WRITE(2,*) 'DATAPACKING=BLOCK'
          WRITE(2,*) ((X(I,J),I=-1,NI+2),J=-1,NJ+2)
          WRITE(2,*) ((Y(I,J),I=-1,NI+2),J=-1,NJ+2)
        CLOSE(2) 
      
      RETURN  
      END SUBROUTINE READER